close all
NumOfStocks = 20;
NumOfBlocks = 2;
%[CovBlocks, Combinations, SampleSizes, Sigma, StockReturns]=GetCovBlocks(0);
load data_2_blocks.mat

% norms = {'fro', 1, 2};
alpha = 0.01;
noise = [4 1];
use_cvx = 0;

CovBlocksNoisy = addNoise(noise, CovBlocks);

iter = 1;
w = cell(1);
w{iter} = (1/length(CovBlocks)) * ones(length(CovBlocks),1); % initial w is feasible
while (true),
    cvx_begin
%        cvx_quiet(true); 
        variable Sigma_hat(NumOfStocks, NumOfStocks) symmetric;
    
        %Define objective function
        f = 0;
        s_max = 0;
        for q=1:NumOfBlocks
            %Add the norm of differences between block covariance and the
            %corresponding block in X.
            f = f + w{iter}(q)*norm(Sigma_hat(Combinations{q}, Combinations{q})-CovBlocksNoisy{q},'fro');
           
            s_max = max(s_max, max(max(CovBlocksNoisy{q})));
        end
        minimize (f)
    
        subject to
        Sigma_hat == semidefinite(NumOfStocks)
        Sigma_hat(:) <= s_max;
        Sigma_hat(:) >= -s_max;
    cvx_end
    
    fs(iter) = f;
    Sigmas_hat{iter} = Sigma_hat;
    error(iter) = norm(Sigma_hat-Sigma, 'fro');
    
    g = [];    
    for i=1:NumOfBlocks,
        g = [g ; norm(Sigma_hat(Combinations{i}, Combinations{i})-CovBlocks{i},'fro')];
    end
    
    x = w{iter} - alpha*g;
    
    if (use_cvx == 0),
        x_sorted = sort(x);
    
        i=[NumOfBlocks-1:-1:0]';
        nu = 1./(NumOfBlocks-i) .* (cumsum(x_sorted(end:-1:1))-1);
        nu = nu(end:-1:1);
    
        L_max = -100;
        for i=1:NumOfBlocks,
            L = 0.5*norm(pos(x-ones(NumOfBlocks,1)*nu(i))-x)^2 + nu(i)*(ones(1,NumOfBlocks)*pos(x-ones(NumOfBlocks,1)*nu(i))-1);
            if (L > L_max )
                L_max = L;
                nu_opt = nu(i);
            elseif (L > L_max & i==NumOfBlocks & x_sorted(i) <= nu(i))
                L_max = L;
                nu_opt = nu(i);
            elseif (L > L_max & i==1 & x_sorted(i+1) >= nu(i))
                L_max = L;
                nu_opt = nu(i);
            end
        end
    
        w{iter+1} = pos(x-ones(NumOfBlocks,1)*nu_opt);
    else
        cvx_begin
        cvx_quiet(true);
        variable y(length(x));
        
        minimize(norm(y-x,2));
        ones(1,length(y))*y == 1;
        y >= 0
        cvx_end;
        w{iter+1} = y;
    end
    
    w{iter+1}(find(w{iter+1}<1e-4)) = 0;
        
    if (iter == 8),
        break;
    end
    
    iter = iter+1;
end

[error_min, ind] = min(error);
w_min = w{ind};
Sigma_hat_min = Sigmas_hat{ind};